In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
def explici_euler(f, x0, times, N):
t0, tf = times
t = np.linspace(t0, tf, N)
dt = (tf - t0) / N
x = np.zeros((len(x0), N))
x[:, 0] = x0
for i in range(N - 1):
x[:, i + 1] = x[:, i] + dt * f(x[:, i])
return t, x
param = (.1, 1.5, .5, .5)
def f(x, t=0):
a, b, c, d = param
I, Z = x
dI = -a * I * Z + b * Z * Z
dZ = -c * Z * I + d
return np.array([dI, dZ])
In [3]:
N = 10000
param = (.1, 1.5, .5, .5)
times = (0, 250)
x0 = (1, 5)
t, x = explici_euler(f, x0, times, N)
plt.plot(t, x.T)
plt.title('I(t), Z(t)')
plt.xlabel('t')
plt.legend(['I(t)', 'Z(t)'])
Out[3]:
In [4]:
n_trial = 100 // 10
Z0 = np.random.normal(5, 5 * .15, n_trial)
I0 = 1
simul = np.zeros((2, N, n_trial))
for k, Z0_ in enumerate(Z0):
_, simul[:, :, k] = explici_euler(f, (I0, Z0_), times, N)
In [5]:
plt.plot(t, np.mean(simul, axis=2).T)
plt.title('mean of I(t) and Z(t)')
plt.xlabel('t')
plt.legend(['mean(I(t))', 'mean(Z(t))'])
plt.figure()
plt.plot(t, np.var(simul, axis=2).T)
plt.title('variance of I(t) and Z(t)')
plt.xlabel('t')
plt.legend(['variance(I(t))', 'variance(Z(t))'])
Out[5]:
In [ ]:
In [ ]: